\documentclass[10pt,a4paper]{article} 

\usepackage{ctex} % 中文支持
\usepackage[top=2.5cm, bottom=2.5cm, left=2.5cm, right=2.5cm]{geometry} % 页边距
\usepackage{amsmath, amssymb} % 数学公式与符号
\usepackage{graphicx}
\usepackage{pythonhighlight}
\usepackage{url} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 使表格美观
\usepackage{array}
\newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
\setlength\extrarowheight{3pt}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{titling}
\setlength{\droptitle}{-2cm} % 标题上移

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{blkarray} % 给矩阵的行列加名称

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%文档的题目、作者与日期
\author{五六七 }
\title{马尔可夫链的几个应用 }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\maketitle

\begin{abstract}
使用马尔可夫链预测机器出现故障的概率，顾客购买不同品牌的味精的选择，以及翻译一篇用替换方法加密的文件。
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\setcounter{tocdepth}{2}
%\renewcommand\contentsname{目录}
%
%\renewcommand {\baselinestretch} {1.3}\normalsize 
%\tableofcontents 
%\renewcommand {\baselinestretch} {1.0}\normalsize


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{问题描述A：计算机工作状态问题}
研究者每隔15分钟观测一台计算机的运行情况。在24小时内收集到97个数据。
其中1表示正常工作状态，0表示不正常工作状态。数据如下：

\vspace{0.3cm}
\begin{center}

1110010011111110011110111111001111111110001101101

111011011010111101110111101111110011011111100111

\end{center}

\vspace{0.3cm}

使用马尔可夫链模型，计算这台计算机的状态转移概率。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{建立模型A}

设 $\{\xi_1,\xi_2,\cdots,\xi_n, \cdots\}$ 是一个随机变量的序列，取值的状态空间 $E$ 为有限集或可数集。
如果对任意正整数 $m,n$ 以及任意状态 $i,j,i_1,\cdots,i_{n-1}$,  都有
\begin{eqnarray}
P\left\{ \xi_{n+m}=j \mid \xi_n=i, \xi_{n-1}=i_{n-1}, \cdots,\xi_1=i_1 \right\} = P\left\{ \xi_{n+m}=j \mid \xi_n=i \right\},
\end{eqnarray}
那么称这个随机变量序列是一个马尔可夫链。 如果上述右端的概率与 $n$ 无关，那么称这是一个时齐的马尔可夫链。称条件概率
\begin{eqnarray}
p_{ij}(m) := P\left\{ \xi_{n+m}=j \mid \xi_n=i \right\}
\end{eqnarray}
为从状态 $i$ 到状态 $j$ 的 $m$ 步转移概率。

在本题中，设 $X_n$ 是第 $n$ 个时段的计算机状态，状态空间为 $E=\{0,1\}$. 为了估计一步转移概率，需要计算观测到的状态数据字符串中字符串 `00', `01', `10', 和 `11' 出现的次数。
编程计算可得如下频数表。

\begin{table}[ht!]\centering
\caption{工作状态转移的频数} \vspace{0.2cm}
\begin{tabular}{|M{3cm}|M{3cm}|}\hline
相邻两个状态 & 频数 \\ \hline 
00 & 8 \\ \hline 
01 & 18 \\ \hline 
10 & 18 \\ \hline 
11 & 52 \\ \hline 
\end{tabular}
\end{table}

由此可得一步转移概率的估计值。
\begin{eqnarray}
p_{00} &=& P\{ X_{n+1}=0 \mid X_n=0 \} = \frac{8}{8+18} = \frac{4}{13}, \\ 
p_{01} &=& P\{ X_{n+1}=1 \mid X_n=0 \} = \frac{18}{8+18} = \frac{9}{13}, \\ 
p_{10} &=& P\{ X_{n+1}=0 \mid X_n=1 \} = \frac{18}{18+52} = \frac{9}{35}, \\ 
p_{11} &=& P\{ X_{n+1}=1 \mid X_n=1 \} = \frac{52}{18+52} = \frac{26}{35}. 
\end{eqnarray}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{编程计算A}

载入数值计算包 numpy. 
\begin{python}
import numpy as np
\end{python}

打开数据文件，保存为一个字符串。
\begin{python}
with open('data15_5.txt') as f:
    s = f.read().replace('\n','')
\end{python}

定义2乘2的零矩阵，用来保存统计数据。
\begin{python}
a = np.zeros((2,2))  
\end{python}

定义一个查找给定子字符串的函数。返回给定子字符串在长字符串中出现的所有位置。
\begin{python}
mfind = lambda s,c: [x for x in range(len(s)) if s[x:x+2] == c]
\end{python}

调用自定义函数，计算 `00', `01' , `10', `11' 四个字符串出现的次数。
\begin{python}
for i in range(2):
    for j in range(2):
        a[i,j] = len(mfind(s,str(i)+str(j)))
\end{python}

打印四个次数组成的2乘2矩阵，打印矩阵每行的和。
\begin{python}
print('matrix a: ', a); 
print('row sum of a:', a.sum(axis=1))
\end{python}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{回答问题A}

根据观测数据，每隔15分钟，计算机从不正常状态转为正常状态的概率是 $9/13$, 停留在不正常状态的概率是 $4/13$. 
而从正常状态转为不正常状态的概率是 $9/35$, 停留在正常状态的概率是 $26/35$. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage 
\section{例子15.6：问题 }
设一个随机系统状态空间为 $E=\{1,2,3,4\}$, 记录观测系统所出状态如下：

\vspace{0.3cm}

\begin{center}
4  3  2  1  4  3  1  1  2  3  2  1  2  3  4  4  3  3  1  1  
1  3  3  2  1  2  2  2  4  4  2  3  2  3  1  1  2  4  3  1
\end{center}

\vspace{0.3cm}

如果这个系统可以马尔可夫链模型来描述，请估计转移概率 $p_{ij}$. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{例子15.6：解答 }

记 $n_{ij}$ 是从状态 $i$ 到状态 $j$ 的转移次数，记 $n_i=\sum\limits_{j=1}^{4}n_{ij}$ 是从状态 $i$ 出发的次数。
则一步转移概率 $p_{ij}$ 可以由 $\frac{n_{ij}}{n_i}$ 来估计。 

载入数值计算包 numpy，载入数据。
\begin{python}
import numpy as np
s = np.array([4, 3, 2, 1, 4, 3, 1, 1, 2, 3,
2, 1, 2, 3, 4, 4, 3, 3, 1, 1,
1, 3, 3, 2, 1, 2, 2, 2, 4, 4,
2, 3, 2, 3, 1, 1, 2, 4, 3, 1])
\end{python}

计算字符串中相邻两个数字出现的次数，保存在矩阵 a 中。
\begin{python}
a = np.zeros((4,4))
for k in range(len(s)-1):
    a[s[k]-1,s[k+1]-1]+=1
print(a)
\end{python}

求得矩阵为 
\begin{eqnarray}
a= \begin{bmatrix} 
 4 & 4 & 1 & 1 \\ 
 3 & 2 & 4 & 2 \\ 
 4 & 4 & 2 & 1 \\ 
 0 & 1 & 4 & 2 \\
\end{bmatrix}.
\end{eqnarray}

对矩阵 a 按行求和，计算从每个状态出发的次数。
\begin{python}
b=a.sum(axis=1)
print(b)
\end{python}

求得从每个状态出发的次数为
\begin{eqnarray}
b = \begin{bmatrix} 
10 \\ 11 \\ 11 \\  7 \\ 
\end{bmatrix}.
\end{eqnarray}

计算从每个状态出发，到达另一个状态的频率。
\begin{python}
p = np.zeros((4,4))
for i in range(4):
    for j in range(4):
        p[i,j] = a[i,j]/b[i]
print(p)
\end{python}

求得一步转移概率矩阵的估计为
\begin{eqnarray}
p = \begin{bmatrix} 
0.4 & 0.4 & 0.1 & 0.1 \\ 
 0.2727 & 0.1818 & 0.3636 & 0.1818 \\ 
 0.3636 & 0.3636 & 0.1818 & 0.0909 \\ 
 0 &  0.1428 & 0.5714 & 0.2857 \\ 
\end{bmatrix}. 
\end{eqnarray}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage 
\section{例子15.7：问题 }
设市场上供应A, B, C三种不同厂家生产的味精，用 `$\xi_n=1$', `$\xi_n=2$', 和 `$\xi_n=3$' 分别表示顾客第 $n$ 次购买A, B, C 厂的味精。设顾客第 $n+1$ 购买时的选择，第 $n$ 次购买的影响覆盖了以往各次购买的影响。
设顾客第一次购买这三个厂家的味精的概率分别为 $0.2, 0.4, 0.4$. 设顾客这次购买对下次购买的影响如下表。

\begin{enumerate}
\item  求顾客第四次购买各家味精的概率。
\item  经过多次购买，三个厂家的市场占有情况如何？
\end{enumerate}

\begin{table}[ht!]\centering
\caption{状态转移概率 } \vspace{0.2cm}
\begin{tabular}{|M{2cm}|M{2cm}|M{2cm}|M{2cm}|}\hline
 & 这次购买A & 这次购买B & 这次购买C \\ \hline
上次购买A & 0.8 & 0.1 & 0.1 \\ \hline
上次购买B & 0.5 & 0.1 & 0.4 \\ \hline
上次购买C & 0.5 & 0.3 & 0.2 \\ \hline
\end{tabular}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{例子15.7：解答 }

根据一步转移概率的定义，可得顾客第二次购买这三个厂家的味精的概率为 
\begin{eqnarray}
P^{(2)} = P^{(1)} P = \begin{bmatrix}0.2 & 0.4 & 0.4 \end{bmatrix} 
\begin{bmatrix} 0.8 & 0.1 & 0.1 \\ 0.5 & 0.1 & 0.4 \\ 0.5 & 0.3 & 0.2  \end{bmatrix} 
= \begin{bmatrix}0.56 & 0.18 & 0.26  \end{bmatrix} 
\end{eqnarray}

重复这个过程，可得顾客第四次购买各家味精的概率为
\begin{eqnarray}
P^{(4)} = P^{(1)} P^3 = \begin{bmatrix}0.2 & 0.4 & 0.4 \end{bmatrix} 
\begin{bmatrix} 0.8 & 0.1 & 0.1 \\ 0.5 & 0.1 & 0.4 \\ 0.5 & 0.3 & 0.2  \end{bmatrix}^3 
= \begin{bmatrix} 0.7004 & 0.136 & 0.1636  \end{bmatrix} 
\end{eqnarray}

经过多次购买，可以理解为求极限 $\lim\limits_{n\to\infty} P^{(n)}$. 程序计算可得 
\begin{eqnarray}
P^{(20)} = P^{(1)} P^3 = \begin{bmatrix}0.2 & 0.4 & 0.4 \end{bmatrix} 
\begin{bmatrix} 0.8 & 0.1 & 0.1 \\ 0.5 & 0.1 & 0.4 \\ 0.5 & 0.3 & 0.2  \end{bmatrix}^{19} 
= \begin{bmatrix} 0.7142 & 0.1309 & 0.1547  \end{bmatrix} 
\end{eqnarray}
因此三个厂家的市场占有率分别为 71.4\%, 13.1\% 和 15.5\%. 

\begin{python}
import numpy as np
P1 = np.mat([0.2, 0.4, 0.4])
P = np.mat([[0.8, 0.1, 0.1],[0.5, 0.1, 0.4],[0.5, 0.3, 0.2]])
P2 = P1 @ P  
print('P2:', P2)
P4 = P1 @ P ** 3 
print('P4:', P4)
P20 = P1 @ P ** 19 
print('P20:', P20)
P30 = P1 @ P ** 29 
print('P30:', P30)
\end{python}

也可以用不变分布来计算。设不变分布为 $\pi=(\pi_1,\pi_2,\pi_3)$, 则有线性方程组 $\pi P =\pi$, 即 
\begin{eqnarray}
\begin{bmatrix}\pi_1 & \pi_2 & \pi_3 \end{bmatrix}
\begin{bmatrix} 0.8 & 0.1 & 0.1 \\ 0.5 & 0.1 & 0.4 \\ 0.5 & 0.3 & 0.2  \end{bmatrix}
= \begin{bmatrix}\pi_1 & \pi_2 & \pi_3 \end{bmatrix}. 
\end{eqnarray}
由此解得不变分布。因为本题中的转移概率矩阵是正则的，所以这个不变分布也是极限分布。


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\section{问题描述B：明文密文问题}
破译文献 \cite{persi-mcmcr} 的图\ref{tu-1}中的密文。

\begin{figure}[ht!]
\centering
\includegraphics[height=6cm, width=14cm]{stanford-stat-prison-letter.png}
\caption{一篇密文 }
\label{tu-1}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{建立模型B}

\subsection{字符替换加密 }

猜测这是一个简单的字母替换加密法，即每个符号对应一个英文字母或标点符号。于是问题转化为寻找对应
\begin{eqnarray*}
f : \{ \text{密文字符集} \} \to \{ \text{英文字符集} \}.
\end{eqnarray*}
考虑每个字符出现的频率，然后找出对应。密文字符的频率只有这个样本，英文字符的频率参照《战争与和平》等英文著作算出。考虑相邻两个字符出现的频率。具体讲，考虑每个字符之后出现某个字符的频率。

接下来计算一阶转移概率，并且定义可信度函数。
设字符 $x$ 之后出现字符 $y$ 的概率为 $M(x,y)$. 这个矩阵通过英文著作统计得到。
定义解密对应 $f$ 的可信度函数为
\begin{eqnarray*}
Pl(f\,) =   \prod_{i} \, M(f\, (s_i), f\, (s_{i+1}) ), 
\end{eqnarray*}
其中 $s_i$ 取遍密文中的每个符号，最后一个除外。

\subsection{一个最简单的例子}

设‘英文字符集’为 $A=\{a,b,c,d,s\}$. 设‘一部’英文著作如下，

\begin{center}
d a b c b s d a a a b c s c d b a b b s
\end{center}

那么可以统计得到一阶转移概率矩阵为
\begin{eqnarray*}
M=
\begin{blockarray}{cccccc}
& a & b & c & d & s \\
\begin{block}{c(ccccc)}
  a & 2/5 & 3/5 & 0 & 0 & 0 \\
  b & 1/6 & 1/6 & 2/6 & 0 & 2/6 \\
  c & 0 & 1/3 & 0 & 1/3 & 1/3 \\
  d & 2/3 & 1/3 & 0 & 0 & 0 \\
  s & 0 & 0 & 1/2 & 1/2 & 0 \\
\end{block}
\end{blockarray}.
\end{eqnarray*}

设密文字符集为 $C=\{1,2,3,4\}$. 设密文为 {\color{blue}2 4 2 1 3 2}. 设密文字符与明文字符的对应为
\begin{eqnarray*}
f: C\to A, f(1)=b, f(2)=a, f(3)=s, f(4)=c.
\end{eqnarray*}
那么这个解密对应的可信度为
\begin{eqnarray*}
Pl(f) &=& M(f(2),f(4))\cdot M(f(4),f(2))\cdot M(f(2),f(1))\cdot M(f(1),f(3))\cdot M(f(3),f(2)) \\
&=& M(a,c)\cdot M(c,a)\cdot M(a,b)\cdot M(b,s)\cdot M(s,a) \\
&=& 0\cdot 0\cdot 3/5 \cdot 2/6\cdot 0 =0. 
\end{eqnarray*}

随机改变解密对应 $f$ 的两个取值， 得到另一个对应 $f_*$, 
\begin{eqnarray*}
f_*: C\to A, f(1)=a, f(2)=b, f(3)=s, f(4)=c.
\end{eqnarray*}


\subsection{MCMC算法}

可信度函数的值越大，则密文字符和明文字符之间的这个对应越可信。
\begin{enumerate}
\item  随机猜测一个对应，设为 $f$.
\item  计算可信度函数值 $Pl(f)$.
\item  在 $f$ 基础上随机交换两个字符的对应结果，得到一个新的对应，设为 $f_*$. 
\item  计算可信度函数值 $Pl(f_*)$; 如大于 $Pl(f)$ 则接受 $f_*$ 为更新的对应。
\item  如 $Pl(f_*)$ 小于 $Pl(f)$ 则以 $Pl(f_*)/Pl(f)$ 概率接受 $f_*$, 以其它概率保留 $f$. 
\end{enumerate}

\subsection{密文明文问题的一个测试例子}

为测试上述算法是否有效，考虑一段明文，如图\ref{tu-2}所示。然后设一个加密对应，将其变成密文。然后测试上述算法。

\begin{figure}[ht!]
\centering
\includegraphics[height=3cm, width=12cm]{hamlet-to-be-or-not-to-be.png}
\caption{测试例子的明文 }
\label{tu-2}
\end{figure}

如图\ref{tu-3}所示，第100步的对应，密文译出来仍然不可理解。到2000步的时候，密文翻译出来可以理解了。继续迭代，对应不再有所变化。

\begin{figure}[ht!]
\centering
\includegraphics[height=6cm, width=12cm]{hamlet-to-be-iterate.png}
\caption{测试例子的迭代过程 }
\label{tu-3}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\section{编程计算B}
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{回答问题B}

上述 MCMC 用于密文明文问题，翻译结果如图\ref{tu-4}所示。

 \begin{figure}[ht!]
 \centering
 \includegraphics[height=6cm, width=14cm]{stanford-stat-prison-letter-decrypted.png}
\caption{翻译出来的明文 }
\label{tu-4}
 \end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{参考文献 }
\begin{thebibliography}{99}

%\bibitem{dingtongren} 丁同仁、李承治，常微分方程教程，高等教育出版社，2022年3月第三版。
\bibitem{sishoukui-2} 司守奎,孙玺菁. \emph{Python数学建模算法与应用}, 国防工业出版社. 2022年1月第1版. 
%\bibitem{hexiaoqun-ara} 何晓群. \emph{应用回归分析(R语言版)}. 电子工业出版社. 2017年7月第1版. 
%\bibitem{dalgaard} Peter Dalgaard 著, 郝智恒等译. \emph{R语言统计入门}. 人民邮电出版社. 2014年6月第1版. 
\bibitem{zhangbo} 张波, 商豪. \emph{应用随机过程}. 中国人民大学出版社.  2016年6月第4版. 
\bibitem{persi-mcmcr}  Persi Diaconis. \emph{The Markov Chain Monte Carlo Revolution}. Bulletin of The American Mathematical Society, 2009, 46(2):179-205.

\end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage 

\section{马尔可夫链的定义}

设 $X$ 是一个有限集合。一个马尔可夫链是一个矩阵 $K(x,y)$, 其中 $K(x,y) \ge 0$以及  $\sum_y K(x,y) = 1$对任意 $x$ 都成立。 
矩阵 $K$ 的每一行是一个概率测度，从状态 $x$ 出发，以 $K(x,y)$ 的概率到达状态 $y$, 下一步以 $K(y,z)$ 的概率到达状态 $z$.  
马尔可夫链的一条实现路径例如：$X_0 = x$, $X_1 = y$, $X_2 = z$, $\cdots$. 于是有
\begin{eqnarray*}
P(X_1 = y\mid X_0 = x) &=& K(x,y), \\
P(X_1 = y, X_2 = z\mid X_0 = x) &=& K(x,y)K(y,z), \\ 
P(X_2 = z\mid X_0 = x) &=&  \sum {}_y K(x,y)K(y,z). 
\end{eqnarray*}
一步转移概率矩阵的 $n$ 次方 $K^n$ 的 $x,y$ 位置表示了条件概率 $P(X_n = y\mid X_0 = x)$.

马尔可夫链 $K(x,y)$ 的稳定分布 $\pi$ 是指 $\pi(x)\mid x\in X$, 满足条件 $\sum {}_x \pi(x) = 1$ 以及等式 
\begin{eqnarray*}
\sum_x\pi (x)K(x, y) = \pi(y). 
\end{eqnarray*}

从这个等式可知， $\pi$ 是矩阵 $K$ 的特征值 $\lambda=1$ 的左特征向量。
马尔可夫链基本定理是说，在适当条件下，矩阵 $K$ 的幂次收敛到这样一个矩阵，这个矩阵的每一行都相同，而且就是稳定分布。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{马尔可夫链的一个简单例子}

设太平洋某小岛上的天气变化是一个马尔可夫链，设转移概率矩阵如下，求它的不变分布和稳定分布。
\begin{eqnarray*}
K=\begin{blockarray}{ccc}
& Rainy & Sunny  \\
\begin{block}{c(cc)}
  Rainy & 3/4 & 1/4 \\
  Sunny & 1/5 & 4/5 \\
\end{block}
\end{blockarray}.
\end{eqnarray*}

解答：
\begin{eqnarray*}
 \begin{pmatrix} \pi_1 & \pi _2 \end{pmatrix}
\begin{pmatrix} 3/4 & 1/4 \\  1/5 & 4/5  \end{pmatrix}
=\begin{pmatrix} \pi_1 & \pi _2 \end{pmatrix} 
& \longrightarrow &
\begin{pmatrix} \pi_1 & \pi _2 \end{pmatrix}
=\begin{pmatrix} 4/9 & 5/9 \end{pmatrix}, \\
%
 \begin{pmatrix} 3/4 & 1/4 \\  1/5 & 4/5  \end{pmatrix}^n
& \longrightarrow & 
\begin{pmatrix} 4/9 & 5/9 \\  4/9 & 5/9  \end{pmatrix}. 
\end{eqnarray*}

问题：编程计算上一页的马尔可夫链的不变分布和稳定分布。

\begin{python}
import numpy as np
K=np.array([[3/4,1/4],[1/5,4/5]])
ev,P=np.linalg.eig(K.T)
myvector=P[:,1]
mypi=myvector/np.sum(myvector)
print('The invariant distribution is ', np.round(mypi,3))

K1=K
for k in range(15):
    K1=np.dot(K1,K)
print('the stable distribution is ', np.round(K1,3))
\end{python}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{马尔可夫链基本定理}

定理：设 $X$ 是一个有限集合。设 $K(x,y)$ 是 $X$ 上的一个马尔可夫链。设存在 $n_0$ 使得 $K_n(x,y) > 0$ 对任意 $n > n_0$ 度成立，则 $K$ 有唯一的稳定分布 $\pi$, 而且当 $n\to\infty$ 时，成立
\begin{eqnarray*}
K^n(x,y)\to \pi(y) \,\,\text{ 对任意的 }\,\, x,y \in X.
\end{eqnarray*}

在实际应用中，状态空间 $X$ 里的状态很多很多，直接从分布 $\{\pi(x)\mid x\in X\}$ 抽样很困难，而从状态 $x$ 转移到状态 $y$ 比较容易。


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{从马尔可夫链角度来看密文明文问题}

在密文明文问题中，状态空间 $X$ 里的元素是密文字符集到明文字符集的各种对应，
\begin{eqnarray}
X = \left\{ f:\,\, {\color{red}\{ +,//,=,*,... \}} \to {\color{blue}\{a,b,...,z,0,1,2,...,9\}} \right\}. 
\end{eqnarray}

设密文有 $m$ 个不同的字符，明文有 $n$ 个不同的字符，则状态空间的元素个数为
 $\vert X\vert = n(n-1)\cdots (n-m+1)$. 这个状态空间有很多状态。 

这个马尔科夫链的稳定分布如下，%其中的 $z$ 是所有 $\{\pi(f)\mid f\in X\}$ 的和，%其中 $M$ 是英语中两个字符相邻的频率矩阵，
\begin{eqnarray*}
\pi (f) = z^{-1} \prod_i M (f\,(s_i), f\,(s_{i+1})).
\end{eqnarray*}

为解决问题而设计的马尔可夫链为 $J(f,f_*)$, 其中
\begin{eqnarray*}
 J(f,f_*) = \left\{ \begin{array}{ll}
1/(n(n-1)\cdots (n-m+1)), & \text{ 如果 } f,f_* \text{只相差两个字符的值，} \\
0, & \text{ 其它情况。}
\end{array}\right.
\end{eqnarray*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{从给定分布进行抽样，Metropolis 算法}

常数 $z$ 定义如下，注意到直接计算 $z$ 非常困难，
\begin{eqnarray*}
z = \sum_f \prod_i M (f\,(s_i), f\,(s_{i+1})). 
\end{eqnarray*}

问题转化为从稳定分布 $\pi(f)$ 里重复抽样。
因为 $z$ 未知，所以 $p(f)$ 未知，所以直接抽样非常困难。
Metropolis 算法解决了这个问题。

\begin{enumerate}
\item  设 $X$ 是一个有限集，设 $\pi(x)$ 是 $X$ 上的一个分布。只需知道 $\pi(x),\pi(x_*)$ 的相对大小。
\item  设 $J(x,y)$ 是 $X$ 上的一个马尔可夫链，设有 $J(x,y) > 0\Leftrightarrow J(y,x) > 0$. 这个分布 $J$ 与目标分布 $\pi$ 没有关系。
\item  Metropolis 算法将 $J$ 转化成另一个马尔可夫链 $K$, 后者的稳定分布正好是目标分布 $\pi$.  
\item  这个马尔科夫链取为
\begin{eqnarray*}
K(x,y) = \left\{ \begin{array}{ll}
J(x,y), & \text{ if } x\neq y, A(x,y)\ge 1, \\
J(x,y)A(x,y), & \text{ if } x\neq y, A(x,y)< 1, \\
J(x,y)+\sum\limits_{z: A(x,z)<1} J(x,z)(1-A(x,z)), & \text{ if } x=y.
\end{array}\right.
\end{eqnarray*}

\item  其中的接受率为
\begin{eqnarray*}
A(x,y) = \frac{\pi(y)J(y,x)}{\pi(x)J(x,y)}.
\end{eqnarray*}

\item  Metropolis 算法的一个解释：
\begin{enumerate}
\item  从 $x$ 出发，依照概率 $J(x,y)$ 选择 $y$;  
\item  如果 $A(x,y) \ge 1$, 那么跳到状态 $y$; 
\item  如果 $A(x,y) < 1$, 那么按照概率 $A(x,y)$ 跳到状态 $y$; 
\item  如果 $A(x,y) < 1$, 按照剩下的概率 $1-A(x,y)$ 停留在状态 $x$. 
\end{enumerate}

\end{enumerate}

注意：常数 $z$ 在第6步的计算中，在分子分母里抵消了。

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{应用马尔可夫链基本定理}

马尔可夫链 $K(x,y)$ 符合等式
\begin{eqnarray*}
\pi(x)K(x,y) = \pi(y)K(y,x).
\end{eqnarray*}
于是由下述计算可知，$\pi$ 是 $K$ 的不变分布，
\begin{eqnarray*}
\sum_x \pi(x)K(x,y) = \sum_x \pi(y)K(y,x) = \pi(y)\sum_x K(y,x) = \pi(y).
\end{eqnarray*}

由基本定理可得，$K^n$ 将收敛到每一行都是 $\pi$ 的一个矩阵。
马尔科夫链进行很多步之后，状态 $y$ 出现的概率是 $\pi(y)$. 
这样 $\pi(y)$ 取值最大的那个 $y$ 就出现了。


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


